c evaluates the T=T_z=1(nn) energies and the corresponding triplet S=1 
cwavefunction. The values M_L=M_S=\mu=0 are used
	program base
	implicit double complex (a-h,o-q,s-z)
	implicit real*8 (r)
	open(unit=7,file='gpr.res',status='unknown')
c reads the angular momentum and parity
	call potdat
	call readsp
	call readcont
	call dis
	call basesp
	call kinetic
	end program base
cc=======================================================================
	subroutine potdat
cc mxrd= maximun dimension for the coordenate of the s.p. wave function
	parameter(mxrd=1000)
	implicit double complex (a-h,o-q,s-z)
	implicit real*8 (r)
cc rolrtol= accuracy requested in the solution for the w.f.
	common/rolrtol/rolrtol
cc rrxmax= max radius at which u(r) is evaluated
cc rixass= rixass=5. at least in medium and heavy nuclei.
cc rrstep= step at which u(r) is evaluated
	common/c1pot/rrxmax,rixass,rrstep
cc rmatch= is the matching radius
cc mnbo= is related to the first integration points
cc icore= icore=1(2) for neutrons(protons)
	common/c2pot/rmatch,mnbo,icore
cc npar=1(-1) for even(odd) parity
cc at,zt,v_0,r_0,a,v_so are the Woods-Saxon parameters  
	common/c3pot/npar(2),at,zt,v_0(2),r_0,a,v_so(2)
	real*8 at,zt,v_0,r_0,a,v_so
cc nwf= is the numbers of points where u(r) is evaluated
	common/nwf/nwf
cc rolc1= factor which changes from k unit to energy ones	
	common/rolc1/rolc1	
cc pot.dat defines the single-particle hamiltonian
	open(unit=5,file='pot.dat',status='old') 
	read(5,*) rrxmax,rixass,rrstep
	read(5,*) mnbo
	read(5,*) rmatch
	read(5,*) icore
	read(5,*) npar(1),at,zt,v_0(1),r_0,a,v_so(1)
	read(5,*) npar(2),at,zt,v_0(2),r_0,a,v_so(2)
	read(5,*) rolrtol
	nwf=rrxmax/rrstep
	if(nwf.gt.mxrd) stop ' nwf.gt.mxrd, I stop'
	write(7,*)' One-particle hamiltonian parameters read from pot.dat'
	write(7,1800) rrxmax
1800	format(' Maximun radius at which u(r) is evaluated -> ',f7.3)
	write(7,1801) rrstep
1801	format(' Step at which u(r) is evaluated -> ',f7.3)
	write(7,1802) rixass
1802	format(' Imaginary part of the asymptotic distance=',f7.3)
	write(7,1803) mnbo
1803	format(' It is related to the first integration points -> ',i4)
	write(7,1804) rmatch
1804	format(' Matching radius -> ',f7.3)
	write(7,1805) icore
1805	format(' 1(2) for neutrons(protons) -> ',i4)
	write(7,*)' W-S parameter for even s.p.s.: at,zt,v0,r0,a,vso '
	write(7,1806) at,zt,v_0(1),r_0,a,v_so(1)
	write(7,*)' W-S parameter for odd s.p.s.: at,zt,v0,r0,a,vso '
	write(7,1806) at,zt,v_0(2),r_0,a,v_so(2)
1806	format(6f7.3)
	write(7,1807) rolrtol  
1807	format(' Accuracy requested to the s.p. wave function -> ',e10.2/)
cc
        itz=icore	
	if(itz.eq.1) then
	   ap=1.0089d0
	else
	   ap=1.0075d0
	endif
	red2m=at*ap/(at+ap)
	rolc1=0.04783258d0*red2m
	
	rbig=r_0*at**(1./3.)
	write(7,*) 'Nuclear radius -> ', rbig
	write(7,*)
	close(unit=5)
	return
	end
cc=======================================================================
cc Reads the single particle discrete states (s.p.d.s): l,2*j,Re(E),Im(E)
	subroutine readsp
cc ibnd= defines all states initily as bound (-1)
cc nspm= max number of discrete states
	parameter (ibnd=-1,nspm=1000)
	implicit double complex (a-h,o-q,s-z)
	implicit real*8 (r)
cc esp= complex energy of the s.p.d.s
cc iysp= indentifies the kind of s.p.d.s (antibound=-2,bound=-1,res=0,scatt=1)
cc jsp= double total angular momentos of s.p.d.s
cc lsp= orbital angular momentum of s.p.d.s
cc nsp= give the total number of discrete states
	common/sp/esp(nspm),itysp(nspm),jsp(nspm),lsp(nspm),nsp
	open(unit=4,file='sp.dat',status='old')
	nsp=0 !initializes the discrete state number 
	do i=1,nspm !initializes all arrays
	 itysp(i)=ibnd
	 esp(i)=0.d0
	 jsp(i)=0
	 lsp(i)=0
	enddo
100	read(4,*) ity,l,jd,rener,renei !reads the s.p.d.s.
	 if(jd.lt.0) goto 110 !stop reading
	 nsp=nsp+1	
	 if(nsp.gt.nspm)stop ' nsp.gt.nspm'
	 itysp(nsp)=ity
	 esp(nsp)=dcmplx(rener,renei)
	 lsp(nsp)=l
	 jsp(nsp)=jd
	 go to 100
110	continue !it concludes to read the s.p.d.s.
	write(7,*)' Discrete states read from sp.dat'
	do i=1,nsp
	 write(7,901) i,itysp(i),lsp(i),jsp(i),esp(i)
	enddo
 901	format(1x,i3,i6,i6,3x,i4,4x,2f10.3)
	write(7,*)
	close(unit=4)
	return
	end
cc=======================================================================
cc It reads from the contour.dat the scattering states that we want include in
cc the s.p.bases and the asociated contuour
cc Each s.p.s. has its own contour
	subroutine readcont
cc nsegm= maximal number of segments which define the curve
cc ncacm= maximal number of s.p. scattering states
cc ngam= maximal number of points in the integration method
	parameter (nsegm=20,ncacm=30,ngam=1000)
	implicit double complex (a-h,o-q,s-z)
	implicit real*8 (r)
cc rolc1= factor which changes from k unit to energy ones
	common/rolc1/rolc1
cc lcac= orbital angular momentum for the scattering state
cc jcac= double total angular momentum for the scattering state
cc ncac= number of scattering state
	common/cont1/lcac(ncacm),jcac(ncacm),ncac
cc zver= vertixes of the contuour
cc ngaus= number of integration point (gauss number) along the segment
cc nsegc= number of segments in the contour
	common/cont2/zver(nsegm+1,ncacm),ngaus(nsegm,ncacm),nsegc(ncacm)
	
	open(unit=25,file='contour.dat',status='old')
	do i=1,nsegm+1 !we initialize the array
	 do j=1,ncacm
	  zver(i,j)=0.d0
	  lcac(j)=0
	  jcac(j)=0
	  nsegc(j)=0
	 enddo
	enddo
	read(25,*) ie !ie=1(0) the contour is in energy(MeV) (k(1/fm)) plane
	ncac=0
 151	read(25,*) l,j !stars reading the s.p. scattering states
	 if(j.lt.0) goto 158 !If j.lt.0 stops reading 
	 ncac=ncac+1
	 if(ncac.gt.ncacm) stop ' ncac.gt.ncacm. STOP '
	 lcac(ncac)=l
	 jcac(ncac)=j
	 nver=0
 156	 read(25,*) rxei,ryei,nga !begin to read the vertices
	  nver=nver+1
	  ngaus(nver,ncac)=nga
	  if(nga.gt.ngam) stop ' nga.gt.ngam. STOP '
	  if((nver-1).gt.nsegm) stop ' (nver-1).gt.nsegm. STOP'
	  if(ie.eq.1) then
	   zver(nver,ncac)=dcmplx(rxei,ryei) !rxei,ryei are complex energy
	  else
	   ckver=dcmplx(rxei,ryei) !rxei,ryei are complex k
	   zver(nver,ncac)=ckver**2/rolc1
	  endif
	  if(nga.lt.0) goto 157 !if nga.lt.0 reads the next state
	  go to 156
 157	 continue
	 nsegc(ncac)=nver-1
	 goto 151
 158	continue
 	write(7,*)' scatt. states'
	do i2=1,ncac
	 write(7,901) lcac(i2),jcac(i2)
	enddo
	  write(7,*)' Contour (the same for all scatt. states)'
	i2=1
	 do i1=1,nsegc(i2)+1
	  write(7,903) zver(i1,i2),ngaus(i1,i2)
	enddo
 901	format(1x,2i4)
 903	format(1x,2f7.3,i8)
	close(unit=25)
	return
	end
cc=======================================================================
cc Finds out which discrete state are inside the contour
	subroutine dis
	parameter(nsegm=20,nspm=1000,ncacm=30)
	implicit double complex (a-h,o-q,s-z)
	implicit real*8 (r)
	common/sp/esp(nspm),itysp(nspm),jsp(nspm),lsp(nspm),nsp
	common/cont1/lcac(ncacm),jcac(ncacm),ncac
	common/cont2/zver(nsegm+1,ncacm),ngaus(nsegm,ncacm),nsegc(ncacm)
cc edis= complex energy of the discrete states in the base
cc itydis= kind of the discrete state in the base
cc jdis= total angular momentum the state in the base
cc ldis= orbital angular momentum the state in the base
cc ndis= numbers of discrete states in the base
	common/dissp/edis(nspm),itydis(nspm),jdis(nspm),ldis(nspm),ndis       
cc ralfa= slope for each segment
cc deltae= vector difference between two vertixes
cc ibe= have the information about which state is in the base
	dimension ralfa(nsegm,ncacm),deltae(nsegm,ncacm),ibe(nspm)
	do i=1,nspm !initialize the arrays
	 ibe(i)=0
	 edis(i)=0.d0
	 itydis(i)=0
	 jdis(i)=0
	 ldis(i)=0
	enddo
	do i1=1,nsegm !initiate the slope
	 do i2=1,ncacm
	  ralfa(i1,i2)=1.d20
	 enddo
	enddo
cc calculates the slope of each segments in each contour
	do 1000 i2=1,ncac 
	 do 1100 i1=1,nsegc(i2)
	  deltae(i1,i2)=zver(i1+1,i2)-zver(i1,i2)
cc if the real part of "deltae" if very small means that the slope if
cc nearly vertical so we leave the initial value of slope
	  if(abs(dreal(deltae(i1,i2))).lt.1.d-12) goto 420
	  ralfa(i1,i2)=imag(deltae(i1,i2))/dreal(deltae(i1,i2))	!else evaluate the slope
 420	  continue
 1100	 continue
 1000	continue
	ndis=0
	do 200 j=1,nsp !Verify which state is in the contour
	 ity=itysp(j)
	 l=lsp(j)
	 jcont=jsp(j)
cc verify if the state which is in sp.dat is in contour.dat too
	 do 376 ilj=1,ncac
	  if(l.ne.lcac(ilj)) goto 376	!if the state isn't in contour.dat
	  if(jcont.ne.jcac(ilj)) goto 376	!go to the next state
	  i2=ilj !indicates the contour for the l,j state
	  goto 388	!if the state is in both file continue
376	 continue
	 goto 199
388	 continue
	 statetype: select case(ity)
	 case(-1) !bound states
	  ndis=ndis+1
	  ibe(ndis)=j
	 case(0) !resonant states
	  rx=dreal(esp(j))
	  ry=dimag(esp(j))
	  nvlast=nsegc(i2)+1
	  if(rx.lt.dreal(zver(1,i2)).or.rx.gt.dreal(zver(nvlast,i2)))
     +     go to 199	!we pass to next state
	  do 300 is=1,nsegc(i2)	!we look for between which segment it is
	   if(rx.ge.dreal(zver(is,i2)).and.
     + 	   rx.le.dreal(zver(is+1,i2))) then
	    if(abs(ralfa(is,i2)).lt.1.d19) goto 190	
	    if(rx.eq.dreal(zver(is,i2))) then
	     write(*,*)' resonance on the vertical contour.STOP', rx,ry
		 write(7,*)' resonance on the vertical contour.STOP', rx,ry
		 stop
	    endif
190	    continue
cc we verify is the state is over or under the segment
	    rycont=ralfa(is,i2)*
     1		(rx-dreal(zver(is,i2)))+imag(zver(is,i2))
	    if(rycont.eq.ry) then	!we check if it is over the contour
	     write(*,*)' resonance on the contour.STOP', rx,ry
	     write(7,*)' resonance on the contour.STOP', rx,ry
	     stop
	    endif
	    if(abs(ry).gt.abs(rycont)) goto 199	!it means that state is out
	    ndis=ndis+1	!it means the state is inside the contour
	    ibe(ndis)=j
	    goto 199
	   endif
300	  continue
	 case(-2) !antibound states
 	  rx=dreal(esp(j))
	  do 400 is=1,nsegc(i2)	!we look for between which segment it is
	   if(dimag(zver(is,i2)).le.0.d0.and.
	1            dimag(zver(is+1,i2)).le.0.d0) goto 399
	    if(rx.le.dreal(zver(is,i2)).and.rx.ge.dreal(zver(is+1,i2)))
	1      then
	     ndis=ndis+1	!it means the state is inside the contour
	     ibe(ndis)=j
	     goto 199
	    endif
 399	    continue
400	  continue
	 end select statetype
 199	continue
 200	continue
	write(7,*)' INPUTS: Discrete states in the single particle bases'
	do i=1,ndis	!we transfer the states in the contour to this array
	 j=ibe(i)
	 edis(i)=esp(j)
	 itydis(i)=itysp(j)
	 jdis(i)=jsp(j)
	 ldis(i)=lsp(j)
	 write(7,901) i,itydis(i),ldis(i),jdis(i),edis(i)
	enddo
 901	format(1x,i4,2x,i4,2x,2i4,6x,2f7.3)
	write(7,*)
	end
cc=======================================================================
cc It builds the complete single particle base: type,l,2j,e and w.f.
	subroutine basesp
	parameter(nsegm=20,ncacm=30,ngam=1000,nspm=1000)
	parameter(nmax=2000,maxs=250,neqd=nmax,mni=3000,mnia=300)
	parameter(mxrd=1000,maxd=500,nkmax=30,nbase=1000,maxc=nbase)
	parameter(ngauss=32)
cc nbase= maximal number in the single particle bases
	implicit double complex (a-h,o-q,s-z)
	implicit real*8 (r)
	common/cont1/lcac(ncacm),jcac(ncacm),ncac
	common/cont2/zver(nsegm+1,ncacm),ngaus(nsegm,ncacm),nsegc(ncacm)
	common/dissp/edis(nspm),itydis(nspm),jdis(nspm),ldis(nspm),ndis
cc eip= complex energy of the bases 
cc ityip= kind of the state of the bases
cc jip= total angular momentum of the bases
cc lip= angular momentum of the bases
cc nip= number of the state in the bases
	common/ip/eip(nbase),ityip(nbase),jip(nbase),lip(nbase),nip
	common/nwf/nwf
cc040415
cc scawf= scawf is u(r), remember, the wave function is u(r)/r
c	common/acsw/scawf(mxrd)
	dimension scawf(mxrd)
cc040415
cc wfmat= it records all wave functions in the single part. bases
	common/wfmat/ wfmat(nbase,mxrd)
cc=======================================================================
cc ksq= will be used for ksquare cal.
	common/ksq/ rgaussx1(ngauss),rgaussw1(ngauss),
     &           cks(nbase),
     &	         ups(nbase,ngauss),us(nbase,ngauss),v_ws(nbase,ngauss),
     &           cweight(nbase)
cc======================================================================
cc040415
cc igam= 1(0) then the subroutine swf calculates a gamow(scattering) state
c	common/rolo/igam
        common/c1pot/rrxmax,rixass,rrstep
	common/c2pot/rmatch,mnbo,icore
	common/c3pot/npar(2),at,zt,v_0(2),r_0,a,v_so(2)
	real*8 at,zt,v_0,r_0,a,v_so
cc rolrtol= accuracy requested in the solution for the w.f.
	common/rolrtol/rolrtol
cc040415
cc rwei= gauss weight calculated which the numerical recipes subroutine gauleg
cc rabsc= gauss factor calculated which the numerical recipes subroutine gauleg
	dimension rwei(ngam),rabsc(ngam)
cc factg= gauss factor for the scattering wave function integration
	dimension factg(nbase)
	do 1000 i=1,nbase	!we intiated the base array
	 ityip(i)=0
	 lip(i)=0
	 jip(i)=0
	 eip(i)=0.d0
 1000	continue
	nip=0
	
cc for k^2 cal.
	cweight=1.d0
	call gauleg(0.d0,rrxmax,rgaussx1,rgaussw1,ngauss)
	
	write(7,*)' INPUTS: Discrete states recalculated'
	do 2000 i=1,ndis	!arranged the discrete states first
	 nip=nip+1
	 ityip(nip)=itydis(i)
	 lip(nip)=ldis(i)
	 jip(nip)=jdis(i)
	 eip(nip)=edis(i)
cc040503
c	 igam=0	!recalculated the discrete state using subroutine swf
	 igam=ityip(nip)
cc040503
	 eipnip=eip(nip)
	 lipnip=lip(nip)
	 jipnip=jip(nip)
cc040415
c	 call swf(lipnip,jipnip,eipnip)
	 call swf(rrxmax,rixass,rrstep,
     1            rmatch,mnbo,icore,
     2            npar,at,zt,v_0,r_0,a,v_so,
     3            rolrtol,
     4            igam,
     5            lipnip,jipnip,eipnip,
     6            scawf,
     7            cks(nip),
     8	          rgaussx1,us(nip,:),ups(nip,:),ngauss,v_ws(nip,:)
     9	          )
cc040415
	 eip(nip)=eipnip
	 print *,lipnip,jipnip,eipnip
	 write(7,901) nip,ityip(nip),lip(nip),jip(nip),eip(nip)
 901	 format(1x,i4,2x,i4,2x,2i4,6x,2f7.3)


	 write(80,901) nip,ityip(nip),lip(nip),jip(nip),eip(nip)
	 do j=1, nwf	!we transfer the wave function to the wfmat matrix
	 write(80,*) real(scawf(j)),aimag(scawf(j))
       enddo

	 do 2100 j=1, nwf	!we transfer the wave function to the wfmat matrix
	  wfmat(nip,j)=scawf(j)
2100	 continue
cc we re-define the complex energy (the gamow program calculates the wave function
cc with around 4% error
2000	continue
	write(7,*)
	do 4000 i1=1,ncac !includes the scattering states in the bases
	 do 3000 i2=1,nsegc(i1)	!it calculate the gauss coefficients
	  ng=ngaus(i2,i1)
	  rag=0.d0	
	  rbg=1.d0
	  call gauleg(rag,rbg,rabsc,rwei,ng)
	  dz=zver(i2+1,i1)-zver(i2,i1)
	  z0=zver(i2,i1)
	  do 3100 j=1,ng !it calculate the scattering energy
	   nip=nip+1
	   if(nip.gt.nbase) stop ' nip.gt.nbase. STOP '
	   ityip(nip)=1
	   lip(nip)=lcac(i1)
	   jip(nip)=jcac(i1)
	   eip(nip)=z0+dz*rabsc(j)
	   factg(j)=sqrt(dz*rwei(j)) !integration factor
cc040503
c	   igam=1 !calculates the scattering state using subroutine swf
	   igam=ityip(nip)
cc040503
cc remember: w.f. is u(r)/r
cc scawf is u(r)
cc the real part of u(r) is positive at the origin
	   eipnip=eip(nip)
	   lipnip=lip(nip)
	   jipnip=jip(nip)
cc040415
c	   call swf(lipnip,jipnip,eipnip)
	   call swf(rrxmax,rixass,rrstep,
     1              rmatch,mnbo,icore,
     2              npar,at,zt,v_0,r_0,a,v_so,
     3              rolrtol,
     4              igam,
     5              lipnip,jipnip,eipnip,
     6              scawf,
     7	            cks(nip),
     8	            rgaussx1,us(nip,:),ups(nip,:),ngauss,v_ws(nip,:)
     9	            )
cc040415

	   cweight(nip)=factg(j)
cc we pass the w.f. times the factor integration to the w.f. matrix
         write(80,901)nip,ityip(nip),lip(nip),jip(nip),eip(nip)
	   do k=1,nwf	
	    wfmat(nip,k)=scawf(k)*factg(j)
	    write(80,*)real(wfmat(nip,k)),aimag(wfmat(nip,k))
	   enddo



3100	  continue
3000	 continue
4000	continue
	write(80,*) nip
	write(7,*)' pole single particle basis'
	do i=1,nip
	if(ityip(i).ne.1)write(7,901)i,ityip(i),lip(i),jip(i),eip(i)
	enddo

      

	return
	end
cc calculation of k^2 operator in berggren basis
	subroutine kinetic
	implicit double complex (a-h,o-q,s-z)
	implicit real*8 (r)
	parameter(mxrd=1000,maxd=500,nkmax=30,nbase=1000,maxc=nbase)
	parameter(ngauss=32)
	common/c1pot/rrxmax,rixass,rrstep
	common/ip/eip(nbase),ityip(nbase),jip(nbase),lip(nbase),nip
	common/ksq/ rgaussx1(ngauss),rgaussw1(ngauss),
     &              cks(nbase),
     &	           ups(nbase,ngauss),us(nbase,ngauss),v_ws(nbase,ngauss),
     &             cweight(nbase)

	dimension cksquare(nip,nip)

        open(unit=100,file='k2.dat',status='unknown')
        open(unit=102,file='u.dat')

cc count nonzero
	ncount=0
        do i=1,nip
	   do j=i,nip
	      if( (jip(i).eq.jip(j)).and.(lip(i).eq.lip(j)) ) then
		 ncount=ncount+1
	      endif
	   end do
	end do
	         
        write(100,*) ncount


	write(102,*) ngauss
        do num=1,ngauss
	   write(102,*) rgaussx1(num),rgaussw1(num)
	end do
	write(102,*) nip
	
	cksquare=0.d0
	do i=1,nip
	   ck2=cks(i)*cks(i)
	   do j=i,nip
	      if(i.eq.j) then
	      cksquare(i,j)=ck2
	      end if
	      if( (jip(i).eq.jip(j)).and.(lip(i).eq.lip(j)) ) then
		 do num=1,ngauss
		    cksquare(i,j)=cksquare(i,j)-rgaussw1(num)*
     &                            (us(i,num)*us(j,num)*
     &                             v_ws(i,num)*cweight(i)*cweight(j)
     &				   )
		 end do
	         write(100,*) i,j,real(cksquare(i,j)),aimag(cksquare(i,j))
	      endif
	   end do

	   write(102,'(1x,i4,2x,i4,2x,2i4,6x,4(f8.5,2x))')
     &     i,ityip(i),lip(i),jip(i),
     &     real(ck2),aimag(ck2),real(cks(i)),aimag(cks(i))
           do num=1,ngauss
	      ctemp_us=us(i,num)*cweight(i)
              ctemp_ups=ups(i,num)*cweight(i)
              write(102,'(5(E16.9,2x))')
     &                   real(ctemp_us),aimag(ctemp_us),
     &                   real(ctemp_ups),aimag(ctemp_ups),
     &                   real(v_ws(i,num))
           end do
	end do
	return
	end

cc=======================================================================
cc Taken from the Numerical Recipes. Calculates the weigts and points
cc for the gauss integration method. Note: It converts the integration
cc limits between (-1,1)
	SUBROUTINE gauleg(x1,x2,x,w,n)
	INTEGER n
	REAL*8 x1,x2,x(n),w(n)
	DOUBLE PRECISION EPS	!is the relative precision
	PARAMETER (EPS=3.d-14)
	INTEGER i,j,m
	DOUBLE PRECISION p1,p2,p3,pp,xl,xm,z,z1
	m=(n+1)/2 !the roots are symetrics, son it calculates the half
	xm=0.5d0*(x2+x1)
	xl=0.5d0*(x2-x1)
	do 12 i=1,m
	 z=cos(3.141592654d0*(i-.25d0)/(n+.5d0))
 1	 continue
	 p1=1.d0
	 p2=0.d0
	 do 11 j=1,n
          p3=p2
	  p2=p1
	  p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
 11	 continue
	 pp=n*(z*p1-p2)/(z*z-1.d0)
	 z1=z
	 z=z1-p1/pp
	 if(abs(z-z1).gt.EPS)goto 1
	 x(i)=xm-xl*z
	 x(n+1-i)=xm+xl*z
	 w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
	 w(n+1-i)=w(i)
 12	continue
	return
	END
